Hot electron mediated desorption rates calculated from excited state potential energy 

surfaces 
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We present a model for Desorption Induced by (Multiple) Electronic Transitions (DIET/DIMET) 
based on potential energy surfaces calculated with the Delta Self-Consistent Field extension of Den- 
sity Functional Theory. We calculate potential energy surfaces of CO and NO molecules adsorbed 
on various transition metal surfaces, and show that classical nuclear dynamics does not suffice for 
propagation in the excited state. We present a simple Hamiltonian describing the system, with 
parameters obtained from the excited state potential energy surface, and show that this model can 
describe desorption dynamics in both the DIET and DIMET regime, and reproduce the power law 
behavior observed experimentally. We observe that the internal stretch degree of freedom in the 
molecules is crucial for the energy transfer between the hot electrons and the molecule when the 
coupling to the surface is strong. 



PACS numbers: 31.15.xr, 71.15.' 



I. INTRODUCTION 
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The advent of femtosecond lasers has initiated ma- 
jor progress in the study of non-adiabatic surface dy- 
namics on a wide range of systems. Photo-induced des- 
orption had already been observed for a few adsorbate 
systems*!* using low-intensity nanosecond laser pulses, 
but high-intensity femtosecond laser pulses have been 
shown to induce desorption in a large class of adsorbate 
systems^ i a i&L 8 i 9 i 10 i and induce chemical reactions which 
cannot proceed by thermal heating. 11 

The mechanism attributed to these reactions is exci- 
tation of substrate electrons by the laser pulse. A sin- 
gle hot electrons can then interact with an initially un- 
occupied adsorbate resonance thus asserting a force on 
the adsorbate nuclei which may then lead to Desorption 
Induced by Electronic Transitions (DIET). Using fem- 
tosecond lasers it is possible to reach high densities of 
excited electrons resulting in a new dominating mecha- 
nism - Desorption Induced by Multiple Electronic Tran- 
sitions (DIMET)* 2 . where several hot electrons interact 
with the adsorbate. 

A different method to produce hot electron based on 
a Metal-Insulator-Metal (MIM) heterostructure was sug- 
gested by Gadzubi^. With an ideal MIM device it is pos- 
sible to tune hot electrons to any desired resonance of an 
adsorbate system and the approach thereby suggests the 
highly attractive possibility of performing selective chem- 
istry at surfaces. Such devices have been constructed and 
characterized** and comprise a promising candidate for 
future hot electron femto chemistry experiments. 

The theoretical framework to describe the non- 
adiabatic dynamics resulting from a hot electron inter- 
acting with an adsorbate is usually based on the con- 
cept of Potential Energy Surfaces (PES). In the Born- 
Oppenheimer approximation the electrons are assumed 
to remain in their ground state and are thus decoupled 
from the nuclei. This allow one to map out a ground 



state PES for the nuclei by calculating the electronic en- 
ergy for each position of the nuclei. Similarly, when an 
initially unoccupied resonance becomes occupied a new 
excited state PES arises which has its minimum at a dif- 
ferent location than the ground state PES and a force is 
exerted on the adsorbate. Several models have emerged 
to deal with non-adiabatic dynamics at surfaces but they 
are usually limited by the difficulty to obtain reliable ex- 
cited state PESs and most theoretical results are based 
on model potentials^********>*2 

An often used method to treat the extreme DIMET 
regime with many contributing electrons is using an elec- 
tronic friction modeh 21 ' 22 i 23 The hot electrons are then 
assumed to thermalize rapidly and the influence of the 
electrons on the adsorbate are treated statistically using 
an electronic temperature which can be several thousand 
Kelvin. The conceptual picture is that of a hot Fermi dis- 
tribution with a tail partially overlapping an adsorbate 
resonance and thereby exerting a force on the adsorbate. 
However, correct calculation of the temperature depen- 
dent friction still requires knowledge of the excited state 
PES. 

The subject of this paper will be the application of 
two-dimensional excited state PESs to calculate desorp- 
tion probabilities. We will be particularly interested in 
the DIET regime where the hot electron has a known en- 
ergy as relevant for the MIM device and the few-electron 
DIMET regime. Although the friction models have en- 
joyed some succes o 24 ' 25 , there is still a need of a mi- 
croscopic non-statistical model of DIMET to test the 
assumption of thermally equilibrated electrons and to 
bridge the gap to the DIET regime. Furthermore, the 
hot electron femtochemistry relevant to the MIM device 
can certainly not be described using an electronic tem- 
perature since all electrons are tuned to a specific energy. 

We start by summarizing the method of linear expan- 
sion ASCF-DFT— used to calculate the excited state 
PESs and note some qualitative features using CO on 
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Pt(lll) as an example. We then discuss the models used 
to obtain desorption probabilities based on the calcu- 
lated potential energy surfaces. First an adiabatic model 
in which the adsorbate jumps between the ground and 
excited state potentials is presented. A general non- 
adiabatic Newns- Anderson like mode l 27 ' 28 is then intro- 
duced and the connection to potential energy surfaces is 
explained. This model with linear coupling has previ- 
ously been solved 2 ^ and applied to the one-dimensional 
desorption problem with model parameters^. We ex- 
tend these results to a two-dimensional adsorbate and 
obtain the non-adiabatic coupling parameters from cal- 
culated excited state potential energy surfaces. In the 
DIET regime the model will be used to show that for 
small excited state lifetimes the main channel of energy 
transfer is the internal degree of freedom, and we em- 
phasize its importance in desorption dynamics. We com- 
pare the calculated desorption probabilities for CO and 
NO on four transition metal surfaces and note some gen- 
eral features of the desorption dynamics. The scattering 
probabilities obtained in the model are then generalized 
to include adsorbates in any vibrationally excited state 
which allow us to extend the calculations to include a 
substrate temperature and to treat the DIMET regime 
within the model. In the appendix it is shown how to ex- 
pand excited states within the projector augmented wave 
formalism and the results and generalizations of scatter- 
ing amplitude calculations are summarized. 



II. POTENTIAL ENERGY SURFACES 

The potential energy surfaces were obtained using 
the code gpa w 30 i 31 which is a real-space Density Func- 
tional Theory (DFT) code that uses the projector aug- 
mented wave methodi 32 ' 33 In all our calculations we used 
the Revised Perdew-Burke-Ernzerhof (RPBE) exchange- 
correlation functional 3 ^ since this has been designed to 
perform well for molecules adsorbed on surfaces, and has 
been shown to perform better than the original PBE 
functional 3 ^ both for isolated molecules^ and for ad- 
sorbed moleculeshammer. For each metal we set up a 
closed-packed surface consisting of three atomic layers 
with the top layer being relaxed. 10 A of vacuum was 
then introduced above the slab and 0.25 monolayer of 
adsorbate molecules relaxed at either top or hep hollow 
site. We then mapped out two-dimensional ground state 
potential energy surfaces in terms of the internal stretch 
and the center of mass (COM) to surface distance coor- 
dinate using 12 irreducible k-points and a grid spacing of 
0.2 A. 

To find the excited state potential energy surfaces we 
use the method of linear expansion ASCF which we 
have published in a previous work 2 ^ and implemented 
in gpaw. In the previous publication we have tested the 
method against inverse photo-emission spectroscopy, and 
found that it performed well for molecules chemisorbed 
on surfaces^ In each step of the self consistency cycle 



an electron is removed from the Fermi level, the density 
of an excited state is added to the total density, and the 
band energy of this state is added to the total energy. 
To get the band energy right we need to expand the ex- 
cited state on the Kohn-Sham orbitals found in each it- 
eration. The method is thus a generalization of the usual 
ASCF where occupations numbers are changed. Instead 
of changing occupation numbers we occupy an orbital 
which is not an eigenstate of the KS Hamiltonian but a 
superposition of eigenstates in such a way that the state 
is as close as possible to the original molecular state. We 
refer to appendix [A] for details on how to do this within 
the projector augmented wave formalism. The excited 
states used in this paper are the anti-bonding 2ir orbitals 
of NO and CO. 

In the previous publication^ we investigated the influ- 
ence of the interactions between neighboring super cells 
for different super-cell sizes, and found that the size- 
dependency of the excitation energy is consistent with 
an electrostatic dipole-dipole interaction. Already for a 
(2 x 2) surface cell, the interaction energy has become 
small, and furthermore this interaction energy will have 
little influence on the slope of the excited-state PES, and 
thus little influence on the calculated desorption rates. 
For this reason, and to keep the calculations manageable, 
we use a (2 x 2) surface cell. 

As an example we show the two-dimensional excited 
state PES superimposed on a ground state PES in the 
case of CO on Pt(lll) top site in Fig. [TJ The molecules 
adsorb with the molecular axis perpendicular to the sur- 
face with carbon closest to the top site. Due to the sym- 
metry of the 27r orbital and the geometry at the ground 
state minimum we cannot induce forces parallel to the 
surface if the molecule is at the ground state minimum 
when excited. The excited state could have unstable 
extremal points with respect to the degrees of freedom 
parallel to the surface, but the model we apply in this 
work only depend on the degrees of freedom with non- 
vanishing derivatives on the excited state PES and we 
thus assume that the center of mass (COM) and internal 
stretch degrees of freedom should capture the essential 
desorption dynamics of the considered systems. 

Since the excited molecule has an extra electron in an 
anti-bonding orbital the excited molecule is expected to 
have a larger equilibrium bond length and this is also 
what we observe. A popular and conceptually simple way 
of explaining desorption in one-dimensional models of 
DIET is the Antoniewicz mechanism^ 7 - where the excited 
molecule induces an image charge on the surface which 
results in an attractive force on the surface. The ex- 
cited molecule is then accelerated towards the surface and 
eventually decays to the steep wall of the ground state 
Morse potential. From Fig. [T]we observe a qualitatively 
different behavior: the COM of the excited molecules ex- 
perience a repulsive force accelerating the COM of the 
molecule away from the surface. This is due to the effect 
of the bond length expansion and the fact that the 2tt 
orbital have a large density in the vicinity of the carbon 
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FIG. 1: Ground and excited state (2tt) potential energy sur- 
faces for CO adsorbed on Pt(lll) - top site. The coverage is 
0.25 monolayer. 



atom which gives a repulsion that dominates the image 
charge attraction. It will be shown below that for the 
considered systems it is primarily excitation of the inter- 
nal degree of freedom which is responsible for the large 
energy transfers leading to desorption. 

The potential energy surfaces for CO adsorbed on Pd, 
Rh and Ru show very similar qualitative features. 



III. MODELS 

The timescale at which adsorbates dissipate energy to 
the substrate is typically on the order of picoseconds 38 
and since the oscillation times for the two modes is 
~ 15 — 50 fs (see section IIV Ap we will assume that 
the molecule has plenty of time to desorb if it has ab- 
sorbed the required energy from a hot electron. This is 
the major assumption we will impose and thus when we 
refer to desorption rates in the following it is the rate of 
transferring at least of the energy needed for a molecule 
to desorb. 

Assuming a Lorentzian resonance with Full Width at 
Half Maximum (FWHM) T = K/t centered at e a , the 
probability that a hot electron of energy e desorbs the 
molecule becomes 



Des 



(e) = 



(r/2) 2 



i 



( £ - £a ) 2 + (r/2)2 



P(t d )e-^ T dt d 



where P{td) is the probability of a desorption event when 
the molecule is excited at t — and decays at time t — 
t d . Using classical dynamics the probability P(td) can 
be obtained by propagating the molecule on the excited 
state PES according to the forces, evaluate the energy 



gain AE after time td and put P{td) = 1 if AE > Eoes 
and P(td) — if AE < Eu es . However, the short lifetime 
of the excited electron implies that classical molecular 
propagation on the excited state PES may not be a good 
approximation. 

In fact, the classical limit is obtained when the action 
S = J dtL(x(t),x(t)) on a representative path satisfies 



\s\ > h. 



(2) 



Assuming a quadratic excited state potential of frequency 
lu and initial potential energy Eq we can evaluate the ac- 
tion on a classical path between initial time U and final 
time tf. For generic timescales one just obtains the usual 
condition of high excitation numbers J?j » Sw, whereas 
for ojAt -c 1, the additional condition of E$At ^> h needs 
to be satisfied in order to apply classical dynamics. In 
the case of CO on Pt(lll) we have E ~ 0.3 eV (Fig. 
[]} and r ~ 1 fs (Fig. |SJ) which gives E At ~ h. Thus 
molecular propagation on the excited state PES is not 
expected to follow the classical equations of motion. Be- 
low we will show an example where a classical analysis 
underestimates desorption probabilities by several orders 
of magnitude. 

This scheme could be extended to a quantum dynami- 
cal treatment of the molecule by propagating the molec- 
ular wavefunction using a two-PES Hamiltonian. How- 
ever, the method still rests on the Born-Oppenheimer 
approximation and the adiabatic concept of potential en- 
ergy surfaces and thus cannot be expected to fully cap- 
ture the non-adiabatic entangled dynamics of the reso- 
nant electron and adsorbate coordinates. 

Instead we consider a Newns- Anderso n 27 ' 28 type 
Hamiltonian with substrate states |fe), a resonant state 
la), adsorbate coordinates adiabatic adsorbate 

ground state potential Vb(xi), and a non-adiabatic cou- 
pling of the resonant electron to adsorbate coordinates 

£a(Xi) 

H = T(±i) + V (xi) + e a {x i )c' a c a 

V " : ^ f V ak ctc k + V: k c{c a ). (3) 



2^ efcCfcCfe 

k 



E 



The strength of the electronic coupling is expressed 
through the function: 



r( £ ) 



k 



\V ak \ 2 6(e-e k ). 



(4) 



The model as such neglects the electron-electron in- 
teraction, but we assume that the important part of 
the electron-electron interactions is the restructuring of 
the metallic electrons when the resonance is occupied 
and that we can capture this effect in an effective non- 
adiabatic coupling. To do this we note that we can ob- 
tain e a {xi) as the expectation value differences of Q 
with the adsorbate at x\ with and without an electron 
in the state \a). Applying this to an interacting problem 
lead us to identify e a (xi) = V\{xi) — Va{xi) where V\{xi) 
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and Vo(xi) are the potential energy surfaces of excited 
and ground states which we have obtained with linear 
expansion ASCF-DFT. 

In the following we will apply the wide band limit 
which means that the individual coupling coefficients V ak 
are assumed to vary slowly in energy and the density 
of states p{e) is taken as constant in the vicinity of the 
resonance. This gives an energy independent coupling 
r = 2np(e a ) J2k \ Vak\ 2 and the resonance spectral func- 
tion corresponding to the electronic part of ([3]) becomes 
a Lorentzian with FWHM T. 

Even in the wide band limit it is quite difficult to han- 
dle the model ([3]) analytically with arbitrary coupling 
function s a (xi). In particular, we would like to calcu- 
late the probability that an incoming substrate electron 
of energy £j scatters inelastically on the resonance and 
is reflected back into the substrate with final energy Ef. 
Fortunately, the potential energy surfaces we are consid- 
ering are close to being quadratic in the region of interest 
(see Fig. [lj and the ground and excited state potentials 
have approximately the same curvature. Taylor expand- 
ing Vb(xi) to second order and e a (xi) to first order in the 
vicinity of the ground state equilibrium positions x® then 
gives 



1. DIET 

In Eq. (|B8|) we show how to calculate the scattering 
probability P ni , nj ( £ i) that an incoming electron of en- 
ergy Si excites the (rij , rij ) mode of a two-dimensional 
harmonic oscillator. The probability of transferring Eji 
or more energy to the adsorbate can then be found by 
calculating 

Pn{si) = ^2 Pn in] {ei)B(hijJini + hw-j-rij - E R ), (7) 
n; , no- 
where 8(x) is the Heaviside step function. The desorption 
rate can then be calculated by integrating this expression 
with the current density of incoming hot electrons. One 
should note that the probability P ni n of exciting the 
(rii,rij) modes in a two- mode model is not just given 
by the product of single mode probabilities P n in a one- 
mode model. This is due to an indirect coupling of the 
two modes through the resonance. The result can be 
generalized to include the substrate temperature and we 
will examine the consequences of this below. 



H = € a c\c a + ^ e k c\c k + (Vakclck + V* k c\c^j 

k k 



(5) 



with e a = Vi{x$) - VbO°) and 



V2dXi 1 



U = 



(6) 



where we have assumed that an appropriate transforma- 
tion to normal coordinates has been performed. Note 
that if the ground and excited state potentials are ex- 
actly quadratic with equal second derivatives, we can re- 
late the coupling constants to the positions x\ of the 
excited state potential minimum as X 2 = hwiAVi with 
AVi = \m l iJ 1 i \x\ — x\\. The quantity = (Xi/hui) 2 
then corresponds to an " initial quantum number" on the 
excited state surface and this becomes the effective di- 
mensionless coupling constant in the model (see appendix 
|B|) . The Hamiltonian (J5]) has previously been subject to 
detailed analysis in the context of inelastic scattering 29 
and applied to desorption dynamics^ for the case of a 
one-dimensional adsorbate with model parameters. 

Below we extend the results of Refs. [H and to a 
two-dimensional adsorbate and calculated the coupling 
parameters Xi from excited state potential energy sur- 
faces. We also calculate scattering amplitudes for an ad- 
sorbate initially in an vibrationally excited state which 
enable us to apply the model to the DIMET regime. 



2. DIMET 

If we assume that the time between individual inelastic 
scattering events is much longer than the scattering time 
itself, it is possible to regard multiple-electron desorption 
events as a sequence of single-electron scattering events. 
Since we have extended the inelastic scattering probabil- 
ities to include situations where the molecule is initially 
in a vibrationally excited state, it is also possible to treat 
DIMET events within the model ([5]). As an example, 
let us assume a single vibrational mode which is initially 
unoccupied (n — 0). When a hot electron with energy 
Ex scatters inelastically on the resonance the result will 
be a probability distribution P 7ll (ei) for all vibrationally 
excited states n of the molecule. If a second electron with 
energy £2 now scatters on the resonance the probability 
distribution will change to P n2 ( e ii £ 2) and so forth. The 
probability P n2 _ ni (£2, ni) of exciting the state tt-2 given 
that the initial state were n\ is calculated in Eq. (|B6[) 
and we can write 



Pn 2 (ei)E2) 



rai=0 



Pn 2 -n 1 (£2,n 1 )P ni (si) 



(8) 



for a two-electron event and similar expressions for 
multiple-electron events. Given an initial distribution 
of hot electrons we may then calculate the probability 
of a desorption event with any number of contributing 
electrons. 
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IV. RESULTS 
A. Parameters 

The parameters in the desorption model ([5]) are the 
width of the resonance T, the frequencies of the normal 
modes Wj, the excitation energy e a and the non-adiabatic 
coupling coefficients A^. We cannot calculate T from first 
principles but we estimate its value from the Kohn-Sham 
projected density of states. It is typically on the order 
of 1 eV, but it will be instructive to treat it as a free 
parameter and examine how it affects desorption proba- 
bilities. 

The frequencies are obtained from a standard normal 
mode analysis and e a is obtained as the excitation energy 
at the ground state potential minimum. The coupling co- 
efficients are determined by mapping out a small area of 
the excited state potential energy surface in the immedi- 
ate vicinity of the ground state potential. In each of the 
considered systems we optimize the area such that it is 
small enough to be linear but large enough to suppress 
numerical fluctuations in the excited state energies. We 
then fit a linear function to this area and transform the 
derivatives to the normal modes. 

In all the considered systems the calculated normal 
modes are similar but not identical to the standard COM 
and internal stretch modes. For example with CO on 
Pt(lll) the internal stretch and COM modes are respec- 
tively: d = (—1,0.75) and z = (1,1) whereas the cal- 
culated modes are in the directions d = (—1,0.68) and 
z = (1,1.11) with respect to the (%c>%o) coordinates 
normal to the surface. Since the desorption probabilities 
are quite sensitive to the value of the non-adiabatic cou- 
pling constants it is important that we take the deriva- 
tives on the excited state PES with respect to the correct 
normal modes. 

Tables U and [IT] below display the calculated param- 
eters. We have only examined CO at on-top sites and 
NO at hep hollow sites. NO is seen to have much lower 
non-adiabatic coupling coefficients and excitation ener- 
gies than CO. The low excitation energies is due to the 
fact that NO already has one electron in the anti-bonding 
orbital and the resonance thus has to lie close to the 
Fermi level of the metal. The small coupling coefficients 
can also be traced to the ground state occupation of the 
2ir orbital on NO. In the Kohn-Sham picture we can 
imagine the resonance corresponding to 2tt lying right 
at the Fermi level being partially occupied. When an ex- 
tra electron is put into the orbital the resonance energy is 
increased due to Hartree repulsion and the initial partial 
occupation is lost. In the true system things are more 
complicated, but the qualitative features are the same: 
exciting NO results in less charge being transferred to 
the molecule than exciting CO and thus a weaker non- 
adiabatic coupling. Thus it is much harder to transfer 
energy to adsorbed NO compared to CO in a one elec- 
tron event, but since the resonance is located much closer 
to the Fermi level a thermal distribution of hot electrons 



Metal 








A z 


A d 


Pt(lll) 


3.89 


0.054 


0.255 


-0.142 


-0.145 


Pd(lll) 


3.64 


0.061 


0.256 


-0.082 


-0.164 


Rh(lll) 


3.80 


0.048 


0.247 


-0.129 


-0.132 


Ru(0001) 


3.74 


0.054 


0.255 


-0.134 


-0.120 



TABLE I: Parameters for CO adsorbed at top site on four 
transition metals. All number are eV. 



Metal 








A z 


A d 


Pt(lll) 


1.71 


0.039 


0.196 


-0.050 


-0.053 


Pd(lll) 


1.48 


0.055 


0.201 


-0.046 


-0.053 


Rh(lll) 


1.82 


0.073 


0.277 


-0.042 


-0.020 


Ru(0001) 


2.14 


0.042 


0.192 


-0.052 


-0.006 



TABLE II: Parameters for NO adsorbed at hep hollow site on 
four transition metals. All number are eV. 



is likely to result in more frequent scattering events than 
for CO. 



B. DIET desorption rates 

The probability that a single electron with energy £j 
scatters inelastically and transfers the energy Er to an 
adsorbate can be calculated in the model ((SJ) with Eq. 
([7]). Our basic assumption is that rate of energy dissipa- 
tion to the substrate is much longer than the time of a 
desorption event and when we refer to desorption rates 
in the following it will mean the rates of transferring the 
energy needed for a molecule to desorb in a truncated 
quadratic potential. 

In Fig. [5] we display the probability that an incoming 
electron will scatter with an energy loss in excess of the 
desorption energy (AE > 1.5eV) for three values of the 
resonance width. When only a single mode is consid- 
ered we see the appearance of oscillator sidebands with 
an energy spacing of hu>. At larger resonance width the 
sidebands are washed out and the probability takes the 
form of a Lorentzian which is detuned by Se a ~ AE/2. A 
simple way to understand this detuning is as a compro- 
mise where both the incoming and outgoing electrons are 
closest to the resonance. Thus we see the emergence of an 
effective inelastic resonance with a center that is detuned 
dependent on the desorption energy and a shape which 
is highly dependent on the lifetime. Such a probabil- 
ity distribution could not have been obtained in a model 
where the transfer of energy to the adsorbate is decou- 
pled from the probability of capturing the electron, and 
the desorption probability would always be a Lorentzian 
(in the wide band limit) centered at e a and multiplied by 
a factor dependent on the details of the potential energy 
surfaces. For T > 0.5 eV the COM degree of freedom 
becomes unimportant and the desorption probabilities 
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Both modes Internal mode 




FIG. 2: Desorption probability for CO adsorbed on Pt(lll) 
for three different values of the resonance width. For V > 0.5 
the one and two-mode probability distributions becomes iden- 
tical Lorentzians with an integrated probability that decays 
exponentially with resonance width (see Fig. |3j). 



obtained using both modes and only the internal degree 
of freedom become identical. 

Assuming an energy independent current of hot elec- 
trons we can integrate the desorption probabilities in 
Fig. [2] to obtain a desorption rate normalized to the 
incident flux of electrons. In Fig. [3] we show how each 
of the two modes contribute to the desorption rate and 
compare with a calculation within the classical adiabatic 
model ([1]). The two single mode rates are obtained by 
setting gd and g z to zero in Eq. (|B8j) . It is seen that it is 
the internal stretch mode that governs the energy trans- 
fer completely in the large width regime and the COM 
mode governs the energy transfer at low width. The rea- 
son for this partitioning is the timescale associated with 
the two different modes. As seen from tables HI and [TTI the 
non-adiabatic coupling constants have approximately the 
same magnitude for the two modes. However, the period 
of oscillation is 5 times larger for the COM mode and for 
small lifetimes there is not enough time to transfer energy 
to the COM mode. From Fig. [3] we see that the maxi- 
mum rate of energy transfer in each mode occurs when 
r ~ hiOi. The desorption rate decrease at small reso- 
nance width, since the hot electron then becomes weakly 
coupled to the resonant state 

In Fig. [5]and [5] we show a comparison of CO and NO 
adsorbed on the different transition metals. Again com- 
paring with tables HI andlTTlit is seen that it is the coupling 
to the internal mode alone which controls the magnitude 
of the desorption rate at large resonance width. Since 
the internal degree of freedom seems to control the rate 
of energy transfer in the physical range of the resonance 
width (typically 0.5 < T < 1.5) we will ignore the COM 
degree of freedom in the following. 



FIG. 3: Desorption rate for CO adsorbed on Pt(lll) as a 
function of resonance width F. In the wide resonance (short 
lifetime) regime the rate is seen to be completely governed by 
the internal stretch excitation whereas the COM excitation is 
governing the desorption rate in the narrow resonance (long 
lifetime) regime. The classical rate becomes several orders of 
magnitude smaller than the quantum rate at large resonance 
width. The inset shows the same data on a logarithmic scale. 
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FIG. 4: Rates of transferring 1.5 eV to CO on 4 transition 
metals. 



1. Comparison of CO and NO 

So far we have analyzed some general features of des- 
orption probabilities and their dependence on the non- 
adiabatic coupling parameters and the lifetime t = h/T. 
Now we will compare the desorption probabilities of CO 
and NO on four transition metal surfaces using experi- 
mentally determined desorption energies. Although sub- 
stantial experimental data exist for various systems in- 
cluding CO and NO, a direct comparison to experimen- 
tal data is difficult since experimental desorption yields 
are highly dependent on the distribution of hot electrons 
in the substrate which depends on the detailed physical 
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pt 

Pd 
Rh 
Ru 




Metal 


E D 


r 


Se 


r>Max 
*D 


Pt(lll) 


1.29 b 


0.8 


0.6 


3 • 10" 11 


Pd(lll) 


LIT 


0.6 


0.6 


5 ■ lO" 9 


Rh(lll) 


1.68 c 


0.4 


0.8 


2 ■ 10~ 15 


Ru(OOOl) 


1.49 d 


0.3 


0.7 


3 ■ 10~ 22 



TABLE IV: Desorption energies and calculated maximum 
desorption probability for NO adsorbed at hep hollow site 
on four transition metals. All numbers except p^ ax is eV. 
'Croci et al^. c Vang et a&>, d Butler et a^ 2 -. 
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FIG. 5: Rates of transferring 1.0 eV to NO on 4 transition 
metals. 



Metal 


E D 


r 


Se 


ryMax 
r D 


Pt(lll) 


1.37" 


1.0 


0.6 


2 ■ 10" 5 


Pd(lll) 


1.48" 


1.5 


0.7 


7- 10~ 7 


Rh(lll) 


1.45" 


1.2 


0.7 


1 ■ 10~ 6 


Ru(0001) 


1.49 a 


0.9 


0.7 


2 • 10 -6 



TABLE III: Desorption energies and calculated maximum 
desorption probability for CO adsorbed at top site on 
four transition metals. All numbers except p^ Iax \ s 
eV. "Experimental values taken from Abild-Pedersen and 
Anderssoo^ 9 ". 



properties of the metal and the applied laser pulse. The 
distribution of hot electrons resulting from a given laser 
pulse could in principle be calculated from first principle, 
however, we will make no attempt of such a calculation 
here but simply compare desorption probabilities of sin- 
gle electron events as relevant for the MIM devic e 13 ' 14 . 
In tables Hill and ITVl we summarize the desorption energy 
Ed, the estimated resonance width T, the detuning of 
the energy at which the incoming electron has the max- 
imum probability of transferring the desorption energy 
Se = e\ lax — e a and the maximum desorption proba- 
bility P^ Iax = PD(ef Iax ) for the four transition metals 
(the maximum probability is detuned from e a as shown 
in Fig. [J). The detuning very nicely follows the rule of 
thumb that Se ~ En/2 in accordance with the picture of 
a compromise between the incoming and outgoing elec- 
tron both being as close as possible to the center of the 
resonance e a . 

In general it is easier for a single electron at the right 
energy to mediate a desorption event involving CO than 
with NO from all the considered systems. However, in 
a femtosecond laser pulse experiment the resulting hot 
electron distribution would have much lower occupation 
numbers at the CO resonances than at a typical NO reso- 
nance. For example taking platinum as an example with 



a thermal electron distribution at 5000 K and referring 
to tables U and |TT] we see that the electronic occupation 
numbers at the resonance energy of CO and NO relates 
as f(^No)/f( e No) ~ 150. We should also note that the 
excited state potential energy surfaces for adsorbed NO 
is only quadratic in a small region near the minimum 
and the model Hamiltonian |5]) is thus not expected to 
describe NO as accurately as CO. 

The desorption probabilities are highly dependent on 
the resonance width T which we can only estimate 
roughly from the Kohn-Sham projected density of states. 
In addition, the electronic lifetime of CO on Pt(lll) has 
been shown to be highly dependent on coverage 3 - since 
the 2tt electrons becomes delocalized and quasistation- 
ary at certain coverages. Furthermore, both CO and NO 
are known to form adsorbate structures which is more 
involve d 41 ' 43 than the simple periodic coverage of 0.25 
monolayer considered here and the dependence of non- 
adiabatic coupling coefficients on coverage certainly de- 
serves a study of its own. 

However, from Figs. |4] and [5] we do observe the gen- 
eral trends that NO has a much weaker non-adiabatic 
coupling to the surfaces than CO and that for both CO 
and NO the coupling to Pt and Pd are similar whereas 
the coupling is weaker for Rh and very low for Ru. This 
decrease in non-adiabatic coupling could hint at a simple 
dependence on the number of d-band electrons. Investi- 
gating this will be the subject of future work. 



C. DIMET desorption rates 

To get an idea of desorption probabilities in the 
DIMET regime we will start by examining how an ini- 
tial excitation influences the probability of transferring a 
given number of vibrational quanta. When the oscillator 
is in an excited vibrational state there is also the possi- 
bility of stimulated emission of vibrational quanta where 
the incoming hot electron gains energy by the scattering 
event. 

In Fig. [5] the maximum probability of transferring An 
quanta is shown for a range of initial quantum numbers n. 
We treat n as a continuous variable since in the case of a 
thermal ensemble of states the initial quantum number is 
simply replaced by a Bose distribution. There is a strik- 
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FIG. 6: Maximum probability of transferring An vibrational 
quanta given that the initial state is n with F = 1.0 eV. 



ing increase in the probabilities of transferring energy to 
the oscillator if the oscillator is already excited. For ex- 
ample, the probabilities of exciting — > 3 and 3 — > 6 are 
3T0 -3 and 2-10 -2 respectively although both transitions 
involve the same energy transfer. Thus if we compare the 
one-electron event Po^6 = 6 • 10~ 6 with the product of 
the two probabilities fb->3->6 = 6 • 10~ 5 we get an order 
of magnitude difference and we still need to include the 
other channels for transferring 6 quanta in a two-electron 
event. 

This also implies that the effect of a finite substrate 
temperature is two-fold: The occupation numbers of ex- 
cited vibrational states will be non-zero meaning that less 
energy transfer is needed to desorb the molecule and the 
likelihood of a given energy transfer is increased if the 
molecule is thermally excited. However at room temper- 
ature the probability that the internal mode is in its first 
excited state is on the order of 10~ 5 and we can safely 
neglect the effect of temperature. 

A hallmark of the DIMET regime is the power law 
dependence of the desorption rate on the laser fluence 
R ~ F n where n depends on the particular adsor- 
bate/substrate system considered^ It is by no means 
trivial that the desorption rate should follow a power 
law and calculating the exponent of a particular system 
is a major challenge of any DIMET model. 

It is reasonable to that assume that the laser fluence 
is proportional to the flux of hot electrons hitting the 
molecule, since the desorption rate typically becomes 
linear— for small fluences corresponding to the DIET 
regime. As a simple model for the desorption rate we 
then consider a given flux J of hot electrons at a fixed 
energy hitting the resonance in equally spaced time 
intervals At = 1/J. We assume that each vibrational 
quantum has a fixed lifetime T v u, and that desorption 
occurs immediately if the vibrational energy reaches the 
desorption energy Ed- The probability that one vibra- 
tional quantum survives the time interval At is e _At / T " ii> 



and the probability of decay is (1 — e~ At / Tvib ). The prob- 
ability that the first electron excites the n'th vibrational 
state is then simply the DIET probability 

Qi(n) = P„(e il 0), (9) 

where P n (si,0) is given by Eq. (|B6p . The probability 
of the adsorbate being in the n'th vibrational state after 
the second electron has scattered is 

oo 

Qzjn) = / J p(m)P n - m (si,m), (10) 

m— 

where P n _ m (e,-,m) is the probability of the transition 
in — ► n (Eq. (|B6|) 1 and p{m) is the probability that the 
adsorbate was initially in the state m given by 

p(m) = £ (e- A */ T ^r(l _ e -At/T ett) *-m 

fc=m ^ ' 

x6(E d -hujk). (11) 

Thus we only sum over values of k below the desorption 
energy since states above Ed would previously have been 
desorbed by assumption. Similarly the probability Qj,{n) 
of being in the n'th excited state after the third scattering 
event can be expressed in terms of Q2(n) and so forth. 
The desorption probability of the N'th electron is then 

P N es =J2^(n)e(tkJn-E D ). (12) 

n 

When enough time intervals is included the probabilities 
converge such that P® es = -Fjv-i an d the desorption rate 
is R{J) = J ■ Pfi es with J = l/At. 

In Fig. [7] we show the rate for NO on Pt(lll) with 
r = 0.8 eV. The desorption energy corresponds to 8 
vibrational quanta. Note that changing the lifetime T v n, 
in this model just corresponds to rescaling the flux. The 
similarity to similar experimental figures^ is striking. At 
small flux the rate is linear whereas it obeys a power law 
(R ~ J n with n > 1) at higher fluences. The fit to 
a power law is very good for fluxes above 0.2 T~}. For 
small values of the detuning (—0.4 < Se < 0.2 eV) we find 
that 5.5 < n < 6 in good agreement with He 44 . For large 
positive values of the detuning the exponent decreases 
dramatically which is probably due the fact that fewer 
transitions dominate the dynamics in this region. This 
means that even though the results was obtained using 
the simple electron flux J(£j) = Jo5(ei—e a —5e) we would 
most likely obtain the same exponent if we generalized 
the model to any flux localized within ±0.2 eV of the 
resonance. 

Although the correspondence with the experimentally 
found exponent may be fortuitous in such a simple model 
the power law itself is very robust to changes in the pa- 
rameters and we obtain similar power laws for CO on 
Pt(lll). For example, changing the value of T results 
in an overall shift of the rates but the exponents are es- 
sentially unchanged. Indeed the exponents appear to be 
determined mainly by the number of vibrational quanta 
needed for desorption. 
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FIG. 7: Desorption rate as a function of electron flux per 
adsorption site. For small electron flux the rate is linear in 
the flux corresponding to the DIET regime whereas for larger 
electron flux the rate obeys a power law (R ~ J" with n > 1) 
corresponding to the DIMET regime. In this figure we show 
the desorption rate of NO on Pt(lll) using the parameters 
given in table ITVl and 7 different values of detuning. 



V. SUMMARY AND DISCUSSION 

We have previously presented a method to obtain 
excited state potential energy surfaces for molecules 
chemisorbed at metal surfaces. 26 In this paper the 
method have been applied and combined with a non- 
adiabatic quantum model to obtain desorption probabil- 
ities for CO and NO on four transition metal surfaces. 

The model we have applied allow us to predict the 
probability that a hot electron will transfer a given 
amount of energy to the different vibrational modes of 
an adsorbate. Our main conclusion is the significant role 
of the internal degree of freedom and the failure of clas- 
sical mechanics to describe the excited state adsorbate 
propagation. 

Combining the model with a simple picture of the de- 
cay and re-excitation of vibrational states reproduces 
the characteristic power laws of DIMET experiments 
and yields the exponent associated with a given adsor- 
bate/substrate system. 

The model we have used for calculating the energy 
transfer rates obviously represents a very simplified view 
of the dynamics. First of all it is a model of non- 
interacting electrons. We assume that we can include 
the important part of the electron-electron interactions 
by using non-adiabatic coupling coefficients Xi obtained 
from the interacting density with linear expansion ASCF- 
DFT. The approximation amounts to assuming ballistic 
hot electrons and instantaneous restructuring of the elec- 
tronic environment when occupying the resonance. Al- 
though this may be the case in some metallic systems 
electron-electron interactions could have effects which go 
beyond a simple renormalization of the non-adiabatic 



coupling. The linear non-adiabatic coupling regime lead- 
ing to Eq. ([5]) corresponds to an assumption of equal 
curvature on the ground and excited state PES. This is 
a good approximation for CO but NO has a very shallow 
excited state PES on some of the transition metals and 
there the approximation may not be as good. 

Furthermore the model assumes that the ground state 
potential is quadratic and that the excited state potential 
is simply a shifted ground state potential. At least in the 
COM direction it is clear From Fig. [1] that the ground 
state potential deviates significantly from a quadratic po- 
tential and since we are concerned with high lying vibra- 
tional excitations this deviation could perhaps have an 
important effect. It may be possible to include anhar- 
monic terms in the Hamiltonian and calculate new scat- 
tering amplitudes perturbatively but the this will be left 
for future work. 

We have focused on the molecules CO and NO, since 
they have a conceptually simple structure and a vast 
amount of experiments have been performed on these 
systems. However, it is well known that GGA-DFT cal- 
culations of CO adsorbed an Pt(lll) predicts CO to bind 
at a hollow site in contradiction to the experimentally 
observed top site 45 . While the difference in adsorbtion 
energy appears to be less with gpaw than in the calcu- 
lations presented in Ref. HH, possibly due to the use of 
the PAW method instead of ultrasoft pseudopotentials, 
the difference is still 80 meV and the inability to pre- 
dict the correct binding site is worrying. On the other 
hand, the existence in the calculation of another adsorp- 
tion site with a slightly lower energy is unlikely to change 
the local shape of the potential energy surface enough to 
qualitatively change the results obtained here. In addi- 
tion, we see a very similar behavior for CO on Ru(0001), 
where DFT does predict the right adsorption site (the 
top site). We have thus chosen to put CO at the ex- 
perimentally observed top site as the hollow site would 
lead to a smaller surface molecule distance and thus very 
different screening and desorption rate. 

As previously mentioned the value of T is estimated 
from the Kohn-Sham projected density of states, but we 
do not know how well this estimate matches the true 
value and as such we have mostly treated T as a free 
parameter. In fact the object of interest in the problem 
is the spectral function of the resonant state, but even 
if we had a reliable way of determining this function we 
would have to make the wide band approximation (where 
the spectral function is a Lorentzian of width T) in order 
to calculate scattering rates. Nevertheless it would be 
very interesting to calculate this function to get an idea 
of the validity of the wide band approximation and to 
obtain a trustworthy value of T. 

We have not made any attempt to predict how the dis- 
tribution of energy evolves after a molecule returns to its 
electronic ground state, but assume that the dissipation 
of energy is slow enough that the adsorbate will desorb 
if the desorption energy has been transferred. This is of 
course a rather crude assumption and the rate of energy 
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transfer should be accompanied by a detailed molecular 
propagation on the full dimensional ground state PES 
to improve the results. Ground state molecular dynam- 
ics would also be necessary to obtain branching ratios 
when there is a possibility of different chemical reactions 
induced by hot electrons. 

However the model we have presented captures some 
of the essential features of non-adiabatic dynamics. For 
example the appearance of an effective inelastic reso- 
nance which is detuned from the electronic resonance by 
an amount depending on the energy transfer is a pure 
non-adiabatic result and would never have emerged from 
an adiabatic model. Furthermore the exponents in the 
DIMET power laws appear to be determined by the num- 
ber of vibrational quanta needed for desorption and thus 
communicates the quantum nature of the dynamics. 

APPENDIX A: PROJECTING KS STATES ON A 
MOLECULAR ORBITAL IN PAW 

The Projector Augmented Wave (PAW) method 32 uti- 
lizes that one can transform single-particle wavefunctions 
\ipn) oscillating wildly near the atom core (all-electron 
wavefunctions), into smooth well-behaved wavefunctions 
\ipn) (pseudo wavefunctions) which are identical to the 
all-electron wavefunctions outside some augmentation 
sphere. The idea is to expand the pseudo wavefunction 
inside the augmentation sphere on a basis of smooth con- 
tinuations \<f>f} of partial waves \(j)f) centered on atom a. 
The transformation is 

llM = \tn) + E (l#> - l#>) (P?l^>. ( A1 ) 
i,a 

where the projector functions \pf) inside the augmenta- 
tion sphere a fulfills 

Y,(Pi\&) = L (m a j) = Sv, \r- R Q | < c 

i 

The method of linear expansion ASCF involves ex- 
panding a molecular orbital \ifi) in Kohn-Sham states 
\ip n ) an d do a self consistent calculation with an addi- 
tional density corresponding to the orbital. 26 . The sim- 
plest way of getting the expansion coefficients is using the 
projector overlaps {ip n \ipi) ~ (ipn\Pi) which is calculated 
in each iteration anyway. However, this method turns 
out to be too inaccurate in the case of CO on Pt(lll) 
due to non-vanishing projector overlaps for highly ener- 
getic Kohn-Sham states as shown in Fig. [8] This implies 
that the expansion coefficients depend on the number of 
unoccupied bands included in the calculation 

To calculate the overlaps (ip n \(pi) exactly, one should 
start by performing a gas-phase calculation of the 
molecule or atom which is to be used in the ASCF cal- 
culation. The pseudo wavefunction ipi(x) corresponding 



to the orbital to be occupied is then saved along with the 
projector overlaps (pt\^i) and the ASCF calculation is 
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FIG. 8: Ground state calculation of CO adsorbed on Pt(lll) 
top site. The projected density of states of the 2tt orbitals 
using the methods of projector/pseudo wavefunction overlap 
and all-electron wavefunction overlap are compared. In the 
projector overlap method the orbital is defined by \p~2-n) = 
■^g(3\px)c — 2\p x )o) which is the orbital most similar to 
the gas-phase calculation. The long high energy tail of the 
projector overlap signals an inaccuracy of the method and 
makes excited state calculations dependent on the number of 
unoccupied bands. Thus we use the all-electron overlaps to 
determine expansion coefficients in this work. 

initialized. In each step of the calculation we can then 
do a numerical integration to obtain the expansion coef- 
ficients by 

Cni = {ipn\il>i) (A2) 

a,j,k 

where (|A1|) was used. Note that there is only a single sum 
over atoms (and only the ones in the molecule) and that 
the cross terms of pseudo/all-electron wavefunction does 
not contribute. This can be seen using the arguments 
following Eq. 20 in Ref. [H 

APPENDIX B: CALCULATING THE INELASTIC 
SCATTERING PROBABILITY 

Here we briefly summarize the calculation leading to 
the inelastic scattering probabilities in the model (JH]). 29 
An explicit expression for the probability has previously 
been obtained 2 ^ for a single mode at initially in the 
ground state. Here we will extend the result to an ex- 
plicit expression for any number of modes initially in a 
thermal ensemble of vibrationally excited states. 

From the Hamiltonian ^ the differential reflection 
matrix R(si,£f) which is defined as the probability per 
unit final state energy that an incoming hot electron with 
energy scatters on the resonance into a final state of £/, 
can be expressed in terms of a four point Green function. 
The inelastic part is contained in the expression: 
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R m (e t ,e f ) = T(e f )T(e l ) 
where the Green functions is 



drdsdt 
2irh 3 ' 



,i[{ ei -e t )r+e t t-e iS yn G ^ g ^ 



G(t, s,t) = 6(s)6(t)(c a (T - S ) C t(T) Ca (i)4(0)>, 



=(*) 



_ AHt/h 



c(0)e" 



iHt/h 



(Bl) 



(B2) 



and () denotes a thermal ensemble of oscillator states. The expression is valid for any non-adiabatic coupling function 
e a (x), but in general it can be very hard to obtain an expression for the Green function. An exception is the wide 
band limit with linear coupling corresponding to the Hamiltonian ((5]). 29 The Green function then becomes 



G(r, s, t) = e{t)9( s )e-^ t - s ^ h - T ^' 2h exp (^gt [i(t - a)ui - (1 + u)fi ~ ruf*] ) 



(B3) 



where e a is center of the resonance, n, is the Bose distribution, g^ = (Xi/htOi) 2 is the effective coupling constant of 
the mode i, and 



fi{r,s,t) 



2 _ e -iu>it _ joiiS + e -iUiT 



(l-e^^l-e™' 8 ). 



(B4) 



The integrals in the scattering matrix (|B1[) can be evaluated by writing the exponentials in (|B3[) as Taylor expansions 
and performing the r integral. This leaves the remaining two integrals as complex conjugates which are evaluated by 
writing factors such as (1 — e luJit ) m by their binomial expansions. For a single oscillator with thermal occupation n 
we obtain the inelastic reflection matrix: 



oo oo 



i? m (e l ,e / ,n)=r 2 e- 2 ^ 1 + 2 ") V V 9 - S(e t - e f - (m, - m 2 )hiu) x F(m, 1 ,m, 2 ), (B5) 

* — ' z — ' mi!m 2 ! 



mi— ni2— 



with 



■F(7nl,m 2 ) 



i=o j=o 



'+:) 



rrii 



7712 

3 



k=Q 1=0 



1 



kin 



(i-j + k-l- g)hw + iT/2 



Although the expression looks rather complicated it has a simple interpretation. Integrating over final state energies 
in the vicinity of An = rni — 777,2 gives the probability of transferring AE = Anhw to the oscillator if the energy of 
the incoming electron is ec 



PAn(si,n) = T 2 e 



-29(1 



i+Wff A "(l + ") 



-F(An, 0) 



.An+l 



(1 + n) An+1 gn 



An! 

5 A "+ 2 (l+n) A "+ 2 ( 5 n) 2 



(An + l)! 
F(An + 2,2) + 



F(An + 1, 1) 



(B6) 



(An + 2)!2! 

where the first term is the probability of adding An bosons, the second term is the probability for removing (coupling 
ng) one and adding (coupling (n + l)g) An + 1 bosons and so forth. 

We can also evaluate the differential reflection matrix for N oscillators initially in the ground state with frequencies 
and coupling constants LOi and gi respectively. The result is 

00 00 mi mw N 

9i 



i? m(£l , £/ ) = r 2 e- 2 £™^ Yl ••■ E 



■9T 



mi— 



mjv— 



mi! . . . m_/v! 



S(sj — Sf — rriihuji 



(B7) 



mi 
h 



m N 
3n 



OO OO l ± 

E-E fl 



■9n 



l 1= l N =0 



hl...l N \ 



£i-e a + iT/2 - J2iLi(3i + h ~ 9i)^ 



It is amusing that the result ()B5[) for a one mode system with initial excitation number n follows from the result (1B7|) 
if we regard (|B5[) as a two-mode system at T = with energies hui and —hoj and coupling constants g(n + 1) and 
gn respectively. For convenience we state the probability of exciting the (m d ,m z ) state from the ground state in the 
two-dimensional model with modes d and z 

n m d m z 

P md mM) = r 2 e- 2 ^+^^-^- (B8) 



m d \n z \ 



EEC- 1 ) 



jd+jz 



3d = ljz 



EE 

k=0 1=0 



9d9i 



1 



fc!Z! e l - e a - (j d + k- g d )^d - {jz + I - 9z)h^ z + iT/2 



12 



1. Elastic scattering 

The elastic part of the scattering matrix for a single oscillator with thermal occupation number n is 

R et (E i ,ef,n)=6(E i -Ef)(l + 2bnG R (e i )), G R (e) = J ^e iEt/h G R {t), (B9) 
G R {t) = -ie(t)(n\c a (t)cl(0)\n). 
We can use the linked cluster theorem to derive the retarded Green function and get the result 



oo oo 



G R {t) = -i0(t) e -9(Wn) e (-ie a -i 9 t^-r/2)t/h y y y n 9 V + n) e -« (ma -m 1 M (BIO) 



TOi!m2! 
mi— m2— 



We can then calculate the elastic part of the scattering probability and get 



PAe- n)-l-r 2 e- g{1+2n) V V 9 9 ( ' - fBll) 

nU^nj-l ^ TOl!m2 , (£ ._ £a _ [m2 _ mi _ 5]M 2 + (r/2) 2- 



mi— 7712— 



When calculating the elastic scattering probability one 
should also remember to include the mi = m,2 terms in 

The n in the expressions above denote the Bose dis- 
tribution and not a specific state |n), but in the context 
of DIMET our main point of interest is the probability 
that a oscillator initially in the state \rii) scatters inelas- 
tically to the state \nj). However, the expression in the 
case of a pure state is very similar to the thermal ensem- 
ble, the only difference being that we should make the 
substitution 

rSW (/ i +/^i nj ( j(/j + / .)) (B12) 

in (|B3p . where L n (x) is the n'th Laguerre polynomial. 
The expression involving Laguerre polynomials is some- 
what more complicated to handle numerically and there- 
fore we have chosen to work with the thermal ensemble 



expressions instead. In the range of parameters in the 
present work the thermal ensemble expressions are very 
good approximations since L n (x) have the same first or- 
der Taylor expansion as e~ nx and for t < h/T we get 
gj, < 0.001. 
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